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Abstract. In |32l 13 1| . fast deterministic algorithms based on spectral methods 
were derived for the Boltzmann collision operator for a class of interactions includ- 
ing the hard spheres model in dimension 3. These algorithms are implemented for 
the solution of the Boltzmann equation in dimension 2 and 3, first for homogeneous 
solutions, then for general non homogeneous solutions. The results are compared 
to explicit solutions, when available, and to Monte-Carlo methods. In particu- 
lar, the computational cost and accuracy are compared to those of Monte-Carlo 
methods as well as to those of previous spectral methods. Finally, for inhomo- 
geneous solutions, we take advantage of the great computational efficiency of the 
method to show an oscillation phenomenon of the entropy functional in the trend 
to equilibrium, which was suggested in the work |16) . 
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1. Introduction 

The construction of approximate methods of solution for the Boltzmann equation 
has a long history tracing back to D. Hilbert, S. Chapmann and D. Enskog ^3] at the 
beginning of the last century The mathematical difficulties related to the Boltzmann 
equation make it extremely difficult, if not impossible, the determination of analytic 
solutions in most physically relevant situations. Only in recent years, starting in the 
70s with the pioneering works by A. Chorin [TJ] and G. Sod j3H], the problem has 
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been tackled numerically with particular care to accuracy and computational cost. 
Even nowadays the deterministic numerical solution of the Boltzmann equation still 
represents a challenge for scientific computing. 

Most of the difficulties are due to the multidimensional structure of the collisional 
integral, since the integration runs on a highly- dimensional unflat manifold. In 
addition the numerical integration requires great care since the collision integral is 
at the basis of the macroscopic properties of the equation. Further difficulties are 
represented by the presence of stiffness, like the case of small mean free path or 
the case of large velocities [TH] . 

For such reasons realistic numerical simulations are based on Monte-Carlo tech- 
niques. The most famous examples are the Direct Simulation Monte-Carlo (DSMC) 
methods by Bird j3j and by Nanbu (35]. These methods guarantee efficiency and 
preservation of the main physical properties. However, avoiding statistical fluctua- 
tions in the results becomes extremely expensive in presence of non-stationary flows 
or close to continuum regimes. 

Among deterministic approximations, perhaps the most popular method is repre- 
sented by the so-called Discrete Velocity Models (DVM) of the Boltzmann equation. 
These methods |2E1 03 El IU3 EE] are based on a cartesian grid in velocity and on a 
discrete collision mechanism on the points of the grid that preserves the main physi- 
cal properties. Unfortunately DVM are not competitive with Monte-Carlo methods 
in terms of computational cost and their accuracy seems to be less than first or- 
der |3b[ H7\ l3*H] . In this work we are interested in high-order deterministic methods 
and therefore we shall not discuss algorithms based on DVM, and we refer the reader 
to the work in preparation |33j . 




More recently a new class of numerical methods based on the use of spectral 
techniques in the velocity space has been developed. The methods were first derived 
in |4L)j . inspired from spectral methods in fluid mechanics jTTj and by previous works 
on the use of Fourier transform techniques for the Boltzmann equation (see j3] for 
instance). The numerical method is based on approximating the distribution func- 
tion by a periodic function in the phase space, and on its representation by Fourier 
series. The resulting Fourier-Galerkin approximation can be evaluated with a com- 
putational cost of 0(n 2 ) (where n is the total number of discretization parameters 
in velocity), which is lower than that of previous deterministic methods (but still 
larger then that of Monte-Carlo methods). 

It was further developed in ^TJ 03] where evolution equations for the Fourier 
modes were explicitly derived and spectral accuracy of the method has been proven. 
Strictly speaking these methods are not conservative, since they preserve mass, 
whereas momentum and energy are approximated with spectral accuracy. This trade 
off between accuracy and conservations seems to be an unavoidable compromise 
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in the development of numerical schemes for the Boltzmann equation (with the 
noticeably exception of [3T?]). 

We recall here that the spectral method has been applied also to non homogeneous 
situations [213 122], to the Landau equation E2], where fast algorithms can be 
readily derived, and to the case of granular gases [33J|2I|. For a recent introduction 
to numerical methods for the Boltzmann equation and related kinetic equations we 
refer the reader to jTSJ . Finally let us mention that A. Bobylev & S. Rjasanow [HIHj 
have also constructed fast algorithms based on a Fourier transform approximation 
of the distribution function. 

In |3*2*1 I3T] a fast spectral method was proposed for a class of particle interac- 
tions including pseudo-Maxwell molecules (in dimension 2) and hard spheres (in 
dimension 3), on the basis of the previous spectral method together with a suitable 
semi-discretization of the collision operator. This method permits to reduce the 
computational cost from 0(n 2 ) to 0(n log 2 n) without loosing the spectral accuracy, 
thus making the method competitive with Monte-Carlo. The principles and basic 
features of this new method will be presented in the next sections. 

The rest of the paper is organized as follows. Section 2 is devoted to a short 
introduction on the Boltzmann equation and its physical properties. Next in Section 
3 we explain the principles of the different spectral algorithms used to compute the 
collision operator. Several numerical results and comparisons to exact solutions 
as well as to Monte-Carlo methods are given in Section 4. An application to a 
challenging non homogeneous test case is finally given in Section 5. Some final 
considerations close the paper in the last Section. 



The Boltzmann equation describes the behavior of a dilute gas of particles when 
the only interactions taken into account are binary elastic collisions. It reads for 
x E f2, v E W 1 where f2 C W 1 is the spatial domain (d > 2) 



where f(t,x,v) is the time-dependent particles distribution function in the phase 
space. The Boltzmann collision operator Q is a quadratic operator local in (t,x). 
The time and position acts only as parameters in Q and therefore will be omitted 
in its description 



2. The Boltzmann equation 



df 
dt 



+ v-V x f = Q(f,f) 



(2.1) 
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In (EH) we used the shorthand / = f(v), /* = f(v*), f = f(v'), f\ = /«). The 
velocities of the colliding pairs (v,v*) and (v',vl) are related by 

1 1 

v' = v — - ((v — v*) — \v — v*\u), vl = v — - ((v — u*) + \v — v*\oS). 

The collision kernel B is a non-negative function which by physical arguments of 
invariance only depends on \v — u*| and cos# = g-u (where g = (v — v*)/\v — In 
this work we are concerned with short-range interaction models. More precisely we 
assume that B is locally integrable. This assumption is satisfied by the hard spheres 
model, which writes in dimension d — 3 

(2.2) B(\v —v*\,cos9) = \v — v*\, 

and is known as Grad's angular cutoff assumption when it is (artificially) extended 
to interactions deriving from a power-law potentials. As an important benchmark 
model for the numerical simulation we therefore introduce the so-called variable hard 
spheres model (VHS), which writes 

(2.3) B(\v — i>*|, cos6>) = C 7 \v — t>*| 7 , 

for some 7 G [0, 1] and a constant C 7 > 0. 

For this class of model, one can split the collision operator as 

Q(fJ) = Q + (fJ)-L(f)f, 

with 

(2.4) Q + (f,f)= [ [ B(\v-v*\,cos9)f'fldadv*, 



(2.5) L(f)= / B(\v -v # \, cos 9) f* da dv # . 

Boltzmann's collision operator has the fundamental properties of conserving mass, 
momentum and energy 



2 



Q(fj) ( j>(v)dv = 0, </>{v) = l,v,\v 

'veR d 

and satisfies well-known Boltzmann's H theorem 

~f flogfdv = -[ Q(fj)\og(f)dv>0. 

Jv£R d Jv£R d 

The functional — J f log / is the entropy of the solution. Boltzmann's H theorem 
implies that any equilibrium distribution function, i.e., any function which is a 
maximum of the entropy, has the form of a locally Maxwellian distribution 

(2.6) M(p,u,T)(v) = _^_exp ( l " & 



(27rT) rf / 2 ^ V 2T 
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where p, u, T are the density, mean velocity and temperature of the gas, defined by 

(2.7) P = [ f( v )dv, u — - I vf(v)dv, T= — f \u — v\ 2 f{v)dv. 

For further details on the physical background and derivation of the Boltzmann 
equation we refer to [13] and [33] . 

3. The spectral methods 

In this section we shall explain the principles of the algorithms to compute the 
collision integral for a fixed value of the spatial variable x. Indeed it is well-known 
that one can reduce to this case by some splitting strategy (see |¥T]l2T)] for example). 

3.1. A general framework. We consider the spatially homogeneous Boltzmann 
equation written in the following general form 

(3-1) % = Q(fJ), 

where Q is given by 

(3.2) Q(fJ)= [ B(y,z)(f'fi-fj)dydz, v G R d 

JUy,z)eC\ 



with 



v' = v + &(y,z), i/ m =v + & tt (y,z), v* = v + Q m (y, z). 

In the equations above, C is some given unbounded domain, and 0, 0', 0^ are 
suitable functions, to be defined later. This general framework emphasizes the 
translation invariance property of the collision operator, which is crucial for the 
spectral methods. We will be more precise in the next paragraphs for some changes 
of variables allowing to reduce the classical operator (j2.1|) to the form ([3.2)1 . 

A problem associated with deterministic methods which use a fixed discretization 
in the velocity domain is that the velocity space is approximated by a finite region. 
Physically the domain for the velocity is M. d , and the property of having compact 
support is not preserved by the collision operator. In general the collision process 
spreads the support by a factor \/2 in the elastic case (see |HJ EOj and also [2H] 
for similar properties in the inelastic case). As a consequence, for the continuous 
equation in time, the function / is immediately positive in the whole domain M. d . 
Thus, at the numerical level, some non physical condition has to be imposed to keep 
the support of the function in velocity uniformly bounded. In order to do this there 
are two main strategies. 

• One can remove the physical binary collisions that will lead outside the 
bounded velocity domain, which means a possible increase of the number of 
local invariants (that is the functions ip such that (ip[ + ip' — <p* — <p) is zero 
everywhere on the domain). If this is done properly (i.e., without removing 
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too many collisions), the scheme remains conservative (and without spuri- 
ous invariants). However, this truncation breaks down the convolution-like 
structure of the collision operator, which requires the translation invariance 
in velocity. Indeed the modified collision kernel depends on v through the 
boundary conditions. This truncation is the starting point of most schemes 
based on Discrete Velocity Models in a bounded domain. 
• One can add some non physical binary collisions by periodizing the function 
and the collision operator. This implies the loss of some local invariants 
(some non physical collisions are added). Thus the scheme is not conservative 
anymore, except for the mass if the periodization is done carefully. In this 
way the structural properties of the collision operator are maintained and 
thus they can be exploited to derive fast algorithms. This periodization is 
the basis of the spectral methods. 

Therefore, we consider the space homogeneous Boltzmann equation in a bounded 
domain in velocity = [T; T] d (0 < T < oo). We need to truncate the integration 
in y and z since periodization would yield infinite result if not. Thus we set y and z 
to belong to some truncated domain Cr C C (the parameter R refers to its size and 
will be defined later). For a compactly supported function with support included in 
Bs, the ball centered at with radius S > 0, one has to prescribe suitable relations 
(depending on the precise change of variable and truncation chosen) between S, R, T 
in order to retain all possible collisions and at the same time prevent intersections of 
the regions where / is different from zero (dealiasing condition). Then the truncated 
collision operator reads 

(3.3) Q R (f, f)= f B(y, z) (£ f - /, /) dy dz 

Jc R 

for v G T>t (the expression for v G M d is deduced by periodization). By making 
some changes of variable on v, one can easily prove for the two choices of variables 
y, z of the next subsections, that for any function ip periodic on T>t the following 
weak form is satisfied 

(3.4) / Q R (f,f)<p(v)dv=- -f [ B(y,z)f.f(<p' m + <p'-<p.-ip)dydzdv. 
Jv T 4 Jv T Jc R 

Now, we use the representation Q R to derive spectral methods. Hereafter, we use 
just one index to denote the d-dimensional sums with respect to the vector k = 
(&4, .., k d ) G Z d , hence we set 

N N 

E== E ■ 

k=-N ki,...,kd=-N 
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The approximate function is represented as the truncated Fourier series 



N 



(3.5) f N (v)= J2 



k=-N 



A = 7^7 / f(v)e-^ k - v dv. 



(2T) £ 



In a Fourier- Galerkin method the fundamental unknowns are the coefficients fk, 
k = —N, . . . , N. We obtain a set of ODEs for the coefficients fk by requiring that 
the residual of ()3.3|) be orthogonal to all trigonometric polynomials of degree less 
than N. Hence for k = —N, . . . , N 

(3.6) jf (^-Q R (f N J N )y-^dv = 0. 

By substituting expression ()3.5|) in ()3.4|) we get 

Q r Un, f N ) = Q R ' + (fN, f N ) - L R {f N ) f N 



with 

(3.7) L R (f N )f N = Yl E P(™>™) fi 



N N 

(l+m)-v 



l=-N m=-N 
N N 

(l+m)-v 



(3.8) Q^+VnJn) = J2 E PMfiLe'* 

l=-N m=-N 

where 

(3.9) (3(1, m)= ! B(y,z)e^( l0 ' {y ' z)+m - @ '' M ) dydz. 



The spectral equation is the projection of the collision equation in IPn, the (2iV + 
l) ci -dimensional vector space of trigonometric polynomials of degree at most N in 
each direction, i.e., 

-j^- = Q R (fN, /at), 

where Vn denotes the orthogonal projection on iFV in L 2 (T>t)- A straightforward 
computation leads to the following set of ordinary differential equations on the 
Fourier coefficients 

df N 

(3-10) ^= E P(l,m)fJ m , 



l-\-m — k 
l,m=-N 
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where (3(1, m) are the so-called kernel modes, given by 

[3{l, m) = [3{l, m) — (3{m, m), 

with the initial condition 

(3.H) A(0) = J^yj v fo(v)e-^dv. 

3.2. Classical spectral methods. In the classical spectral method |H], a simple 
change of variables in (|2.1j) permits to write 

(3.12) Q(f,f)= [ [ B c (g,u)(f(v')f(v'J-f(v)f(v*))dLjdg, 

with g = v — v* G M. d , to G and 

{v' = v - \{g - \g\uS) =: v + &(g, u), 
v'*=v-\{g+\g\u) =:v + &^g,u), 
v* = v + g =: v + B*(g,u). 

Finally B c is defined by 

(3.14) B c (g,iu) = 2 d - 1 (l - (g ■ ^f 2 ' 1 B(\g\, 2(g • uof - l). 

The Boltzmann operator ()3.12j) is now written in the form ()3.2j) with (y, z) = 
(g, oj) G M. d x S^ 1 =: C. Moreover, from the conservation of the momentum v[ +v' = 

+ v and the energy |t>^| 2 + \v'\ 2 = |t>*| 2 + |t>| 2 , we get the following result |4U] . 
assuming supp / C Bs, 

• we have supp Q(f, f) C B^ s , 

• the collision operator is then given by 

Q(fJ)(v)= [ [ B(\g\, cos 9) (f(v')f(vi)-f(v*)f(v))dudg, 

with v', v'^jV* G B (2+V 2 )R . 

As a consequence of this result, in order to write a spectral approximation to 
()3.1|) we consider the distribution function / restricted on [—T,T] d , (0 < T < +oo), 
assuming /(f) = on [—T,T] d \ Bs, and extend it by periodicity to a periodic 
function on [— T, T] d . We truncate the domain for (y, z) = (g, uj) as Cr = Br x § d_1 
for R > (defining thus Q R ). Following the previous discussion on the dealiasing 
condition, we take R = 2S and the shortest period can be restricted to [—T,T] d , 
with T > (3 + y/2)S/2 (see for a more detailed discussion |4*T1 ). 
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Then, we apply the spectral algorithm (J3.7|) and (J3.8|) and get the following kernel 
modes j3 c (l, m) 

(3.15) /3 c (/, m )=/ / B(\g\,cosO)e- i %( 9 -^- i M"-^) dudg. 

Jb r Js*- 1 

We refer to [ITj |2"2] for the explicit computation of Fourier coefficients j3 c (l,m) in 
the VHS case where B is given by (j2.3jl . Now, the evaluation of the right-hand side 
of (|3.1U|) requires exactly 0(N 2d ) operations. We emphasize that the usual cost for 
a DVM method based on N d parameters for / in the velocity space is 0(N 2d M) 
where M is the numbers of angle discretizations. 

3.3. Fast Spectral methods (FSM). Here we shall approximate the collision 
operator starting from a representation which conserves more symmetries of the 
collision operator when one truncates it in a bounded domain. This representation 
was used in |2S] to derive finite differences schemes and it is close to the classical 
Carleman representation (cf. [121 )• The basic identity we shall need is (for u G M, d ) 

(3.16) \! F{\u\a-u)da=-^ ! 5(2y ■ u + \y\ 2 ) F(y) dy. 
2 Jsd-i \u\ d 2 J Rd 

Using ()3.16j) the collision operator (j2.1j) can be written as 

(3.17) Q(fJ)(v)=2 d - 1 [ [ Bf(y,z)5(yz) 

Jx€R d Jy£R d 

(f(v + z)f(v + y)-f(v + y + z)f(v)) dy dz, 

with 

Bf(y, z) = 2 d ^ B f\y + z\, - V ' ( * + Z } ) \v + ■ 

\ \v\\v + AJ 

Thus, the collision operator is now written in the form (|3.2j) with (y, 

C, B(V, Z ) = ^ f (y,z)5(y-z), and < = v + z =: v + & m (y,z), v' = v + y =: v + &(y,z), 

v* = v + y + z =: v + 0*(y, z). 

Now we consider the bounded domain T>t = [T,T] d , (0 < T < oo) for the 
distribution /, and the bounded domain B R x B R for (y, z) (for some R > 0). If / 
has support included in B s , S > 0, geometrical arguments similar to the one for the 
classical spectral methods (see [HJ [321 EI!) show that we can take R = \/2S and T 
as in the classical spectral method to get all collisions and prevent intersections of 
the regions where / is different from zero. The (truncated) operator now reads 
(3.18) 

Q R (f,f)(v)= [ [ B f (y,z)6(yz){f(v + z)f(v + y)-f(v + y + z)f(v))dydz, 

JyeB R JzeB R 

for v G T>t- This representation of the collision kernel yields better decoupling 
properties between the arguments of the operator. From now, we can apply the 
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spectral algorithm (J3.7|) and (J3.8j) to this collision operator and the corresponding 
kernel modes are given by 

p f (l,m)= f [ B(y,z)5(yz)e^( ly+m - z ) dydz. 

JyeB R JzeB R 

In the sequel we shall focus on /3* , and one easily checks that (3^(1, m) depends only 
on \l\, \m\ and \l ■ m\. 

Remark 3.1. Note that the classical spectral method originates the following form 
of the kernel modes in the y, z notation 

PV,m)= [ [ Bf(y,z)5(y z) Xly+zl < R e^( ly+m - z ) dydz. 
JyaB R JzeB R 

One can notice that the condition \y + z\ 2 = \y\ 2 + \z\ 2 < R 2 couples the modulus 
of y and z, such that the ball is not completely covered (for instance, if y and z are 
orthogonal both with modulus R, the condition is not satisfied, since \y + z\ = \JlR). 
This explains the better decoupling properties between the argument of the collision 
operator of this representation. 

4. Fast algorithms 

The search for fast deterministic algorithms for the collision operator, i.e., algo- 
rithms with a cost lower than 0(N 2d+e ) (with typically e = 1 for DVM, or e = 
for the classical spectral method), consists mainly in identifying some convolution 
structure in the operator (see for example [E1E2])- The aim is to approximate each 
(3^(1, m) by a sum 

A 

(4.19) p f (I, m) = J2 a p (l)a' p (m). 

P =i 

This gives a sum of A discrete convolutions and so the algorithm can be computed 
in 0(A iV log 2 N) operations by means of standard FFT techniques To this 
purpose we shall use a further approximated collision operator where the number of 
possible directions of collision is reduced to a finite set. We start from representation 
f!3.18|) and write y and z in spherical coordinates 

Q R (fJ)(v) = ][ [ 6(e-e')dede' 




p d - 2 (p') d - 2 B f (p,p') (f(v + p'e')f(v + pe)-f(v + p'e> + pe)f(v))dpdp> 



(note that thanks to the orthogonality condition imposed by the Dirac mass on y and 
z, 13$ depends only on the modulus of y and z). Let us denote by A a discrete set of 
orthogonal couples of unit vectors (e, e'), which is even, i.e., (e, e') G A implies that 
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(— e, e'), (e, — e') and (— e, — e') belong to ^4 (this property on the set A is required 
to preserve the conservation properties of the operator). Now we define to be 



Q RA (fJ)(v) 

R r R 



dA 



p d ~ 2 (p') d ~ 2 B^p,p') (f( v + p 'e')f(v + pe)-f(v + p'e' + pe)f(v))dpdp' 

-RJ-R 

where dA denotes a discrete measure on A which is also even in the sense that 
dA(e,e') = dA(—e,e') = dA(e,—e') = dA(—e,—e'). It is easy to check that Q R ' A 
has the same conservation properties as Q . We make the decoupling assumption 
that 

(4.20) Vyi_z, B f (y,z)=a(\y\)b(\z\). 

This assumption is obviously satisfied if B* is constant. This is the case of Maxwellian 
molecules in dimension d = 2, and hard spheres in dimension d = 3 (the most rele- 
vant kernel for applications). Extensions to more general interactions are discussed 
in j22]. 

Let us describe the method in dimension d = 3 with satisfying the decoupling 
assumption (|4.2U|) (see |32] for other dimensions). First we change to spherical 
coordinates 



P{l,m) = \ 



5(e ■ e') 



eeS 2 Je'dS 2 



\p\ a(p) dp 



R 



R 



\p'\b(p') e'T'W) dp' 



R 



de de' 



and then we integrate first e' on the intersection of the unit sphere with the plane 

1 



where 



J R,a\ 



65 2 



FrJ} ■ e) 



e'&S 2 r\e A 



^R,b(™ ■ e') de' 



de, 



|p|a(p)e^ s dp, 4>% b (s) 



R 



\p\b(p)e l $ ps dp. 



-R 



Thus we get the following decoupling formula with two degrees of freedom 



p R (l,m) 



0LG-e)V|, 6 (n e x(m))rfe, 



where S 2 , denotes the half-sphere and 



^ (n e x(m)) 



^ b (\U e ,(m)\cos9) d6, 
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and LT e ± is the orthogonal projection on the plane e L . In the particular case where 



Now the function e i— > 4? Ra {l • e) %l)\ b {Yi e ^{m)) is periodic on S% ■ Taking a spherical 
parametrization (9, (p) of e G S\_ and taking for the set A uniform grids of respective 
size Mi and M 2 for 9 and (p we get 



and (6 p ,<p q ) = (p^/M 1 ,qrr/M 2 ). 

We consider this expansion with M = Mi = M 2 to avoid anisotropy in the 
computational grid. By using the Fast Fourier Transform (FFT) algorithm, the 
computational cost of the algorithm is then 0(M 2 N 3 log N), (compared to 0(M 2 N 6 ) 
of a direct discretization on the grid for a DVM method, and 0(N 6 ) of the classical 
spectral method). 

Let us finally mention that the mathematical analysis of the fast algorithm in [32] 
provides the following results: 

• it is spectrally accurate according to the parameters N and M; 

• the error on the conservation laws of momentum and energy is spectrally 
small according to the parameter N, and no additional error (according to 
the speed-up parameter M) is made. 

This two properties were the main motivation for the development of the method 
of described above to obtain the decomposition (j4.19j) . Other advantages are 
that this particular decomposition does not introduce instability in the equation 
(see [321 Theorem 3.1] for instance) and it is naturally adaptative (as it is based on 
the rectangular quadrature rule for approximating integrals of periodic functions). 
Finally, another advantage of the proposed method is that it is still easy to imple- 
ment since it is only based on FFT. 



In this section we will present several numerical results for the space homogeneous 
equation which show the improvement of the fast spectral algorithms with respect 
to the classical spectral methods and how they compare with Monte-Carlo methods. 
The time discretization is performed by suitable Runge-Kutta methods. 





p,q=0 



where 




5. Numerical results in the homogeneous case 



5.1. Spatially homogeneous Maxwell molecules in dimension 2. 
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Comparison to exact solutions. We consider 2D pseudo-Maxwell molecules (i.e., the 
VHS model with 7 = 0). In this case we have an exact solution given by 



with S = 1 — exp(— 1/8)/2, which corresponds to the well known "BKW" solution 
[3]. This test is performed to check spectral accuracy, by comparing the error at a 
given time, when using n v = 8, 16 and 32 Fourier modes for each coordinate. We 
present the results obtained by the classical spectral method and the fast spectral 
method with different numbers of discrete angles. 

Figure ^ shows the relative L°°, L , and L 2 norms of the difference between the 
numerical and the exact solution, as a function of time. These errors are computed 
according to the following formula 



with i = (11,12) and iV = n v /2 for p = 1 and p = 2. A similar expression is 
used for the L°° error. Note that the error increases initially, and then decreases 
almost monotonically in time. After a long time the error starts increasing again. 
This effect is due to aliasing. Indeed, for a fixed computational domain, when the 
number of Fourier modes increases, the effect of aliasing becomes dominant over the 
error due to the spectral approximation. For this reason, the size of the domain 
is chosen in order to minimize the aliasing error. A trade-off should be obtained 
between aliasing and spectral error, which means that the size of the domain should 
be increased when increasing the number of Fourier modes. Roughly speaking, the 
period should be chosen in such a way that the two contributions of the error are 
of the same order of magnitude. In this test, the radius of the ball, which defines 
the computational domain is T = 4 for n v = 8, T = 5 for n v = 16 and T = 7 for 
n v = 32. We refer to [IT] for a more detailed discussion about aliasing. 

Concerning the comparison between the classical and fast spectral methods, we 
observe that for a fixed value of n v , the numerical error of the classical spectral 
method and of the fast algorithm is of the same order. Moreover, the influence of 
the number of discrete angles is very weak. Indeed, with only M = 4, the results 
are quite similar even for large n v and as expected the number of discrete angles 
does not affect the variations of energy, which are of the same order of magnitude as 
the numerical error (note that there is no variation for the momentum since in the 
special case of even solutions, it is preserved to by the spectral scheme). In Table 

we give a quantitative comparison of the numerical error £\ at time T en d = 1. We 
can also observe the spectral accuracy for the classical and fast methods: the order 




exp(-v 2 /2S) 
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Number of 
points 


Classical 
spectral 


Fast spectral 
with M = 4 


Fast spectral 
with M = 6 


Fast spectral 
with M = 8 


8 


0.02013 


0.02778 


0.02129 


0.02112 


16 


0.00204 


0.00329 


0.00238 


0.00224 


32 


1.405E-5 


2.228E-5 


1.861E-5 


1.772E-5 



Table 1. Comparison of the L 1 error in 2D between the classical 
spectral method and the fast spectral method with different num- 
bers of discrete angles and with a second-order Runge-Kutta time 
discretization at time T end = 1. 



of accuracy is about 3 between 8 and 16 grid points, whereas it becomes 7 between 
16 to 32 points. 
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Figure 1. 2D homogeneous case: evolution of the numerical L 1 and 
L°° relative error of f(t,v). 



Efficiency and accuracy. Now, we still consider 2D pseudo-Maxwell molecules (i.e., 
7 = 0) with the following initial datum 

where vq = (1,2). In this case, we do not know the exact solution but we want 
to study the influence of the number of discrete angles on a non-isotropic solution. 



cxp 



\v - v \ 



+ exp 



\v + v \ 



v e 
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Thus, this test is used to check the energy conservation and the evolution of high- 
order moments of the solution. 

The time step is chosen small enough to reduce the influence of the time discretiza- 
tion, i.e., At = 0.025. Moreover, the computational domain is taken large enough 
with respect to the number of grid points in order to reduce the aliasing error due 
to the periodization of the solution. Simulations are performed with n v =16, 32 and 
64 points. 

In Figure El the relaxation of the entropy and the temperature components for 
the fast and classical spectral methods is shown. The energy is conserved by the 
continuous collision operator, but using the spectral method the total energy can 
change with time, it is then a good indicator on the accuracy of the numerical 
solution. Indeed, the total discrete energy is not exactly conserved over time, but 
if aliasing error is small, it is conserved within spectral accuracy, typically here the 
variations are about 10~ 4 when n v = 32. Moreover, we observe that the number of 
discrete angles does not affect too much the transient regime. For instance, with only 
four angles on the half sphere, the relaxation of entropy and temperature components 
are very close to the numerical solution obtained by the classical spectral method. 
Finally, we plot in Figure E3 the time evolution of high-order moments of fw(t, v) 
given in discrete form by 



High-order moments give information on the accuracy of the approximate distribu- 
tion function tail. Once again, we observe that the number of angles does not affect 
the results even if the solution is non-isotropic. 

To conclude, we observe that in dimension d = 2, the fast algorithm is really 
efficient in terms of accuracy and computational cost compared to the classical 
spectral method. In Table 121 we report the computational times of the methods 
which show a speed-up of the fast solver independently of the number of points 
used in our tests and with a maximum speed-up reached for N = 64 where the fast 
methods with M = 4 is more than 17 times faster than the classical method. 

5.2. Spatially homogeneous hard spheres in dimension 3. In this section we 
consider the 3D Hard Sphere molecules (HS) model. The initial condition is chosen 
as the sum of two Gaussians 



with cr = l and Vq = (2, 1,0). The final time of the simulation is T en d = 3 and 
corresponds approximatively to the time for which the steady state of the solution 



N 



M k (t)=Av 2 \vi\ k f N (t,vi). 



l=-N 
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-3.5 
-3.6 
-3.7 
-3.8 
-3.9 
-4 
-4.1 



H with fast spectral T=1 8, M=4, nv=32 + 
H with classic spectral T=18, nv=32 
H with classic spectral T=24, nv=64 
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Figure 2. 2D homogeneous case: relaxation of the entropy and the 
temperature components for the fast and classical spectral methods with 
respect to the number of modes per direction n v and the length box T. 
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Figure 3. 2D homogeneous case: time evolution of the variations of 
high order normalized moments M.^, M.^ and M.§ of f(t,v) for the 
fast and classical spectral methods with respect to the number of modes 
per direction n v and the length box T. 
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Number of 
points 


Classical 
spectral 


Fast spectral 
with M = 4 


Fast spectral 
with M = 6 


Fast spectral 
with M = 8 


16 


2 sec. 40 


1 sec. 15 


1 sec. 70 


2 sec. 30 


32 


38 sec. 01 


5 sec. 55 


8 sec. 47 


11 sec. 10 


64 


616 sec. 


35 sec. 50 


54 sec. 66 


71 sec. 27 



Table 2. Comparison of the computational time in 2D between the 
classical spectral method and the fast spectral method with different 
numbers of discrete angles and with a second order Runge-Kutta time 
discretization. 



Number of 
points 


Classical 
spectral 


Fast spectral 
with M = 4 


Fast spectral 
with M = 6 


Fast spectral 
with M = 8 


16 


1 min. 14sec. 


3 min. 31sec. 


7 min. 45 sec. 


13 min. 44 sec. 


32 


118 min. 02 sec. 


50 min. 31sec. 


105 min. 19 sec. 


186 min. 18sec. 


64 


125/i54 min. 


8/i45 min. 22sec. 


21/i39 min. 


35/i01 min. 28sec. 



Table 3. Comparison of the computational time in 3D between the 
classical spectral method and the fast spectral method with different 
numbers of discrete angles and with a second-order Runge-Kutta time 
discretization. 



is reached. The time step is At = 0.1 and the length box is taken as T = 12 when 
n v = 16 and T = 15 when n v = 32. 

This test is used to check the evolution of moments and particularly the stress 
tensor Pjj, i,j — 1, ■ ■ ■ ,3 defined as 

p i,j = /(«)(«»- -Uj)dv, (i, j) £ {1, 2, 3} 2 , 

where (ui)i are the components of the mean velocity. As in the previous case, 
we compare the classical and fast methods in terms of computational time (see 
Table EJ) and accuracy. In Figure |H we propose the evolution of the temperature 
for the two methods using 32 grid points in each direction. The solution is also 
compared with the solution obtained from the Monte-Carlo method. The discrete 
temperatures agree well in this case and the efficiency of the fast algorithm is verified 
since the computational time is highly reduced using only Mi = M 2 = 4 discrete 
angles without affecting the accuracy of the distribution function. We remark that 
in dimension d = 3 the speed-up of the methods becomes really evident for large 
values of N. Again for iV = 64 and M = 4 the fast methods is more than 14 times 
faster. 
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Figure 4. 3D homogeneous case: comparison between the fast and 
classical spectral methods and the Monte-Carlo methods for the tem- 
perature components relaxation. 



Comparison with Monte-Carlo. Finally we compare the results obtained with our 
spectral method with those obtained by a Monte-Carlo scheme. We use a stan- 
dard version of Monte-Carlo, which may be referred to as the Nanbu-Babowsky 
scheme [US |2] . 

In the case of Monte-Carlo methods, the moments are computed by using unbiased 
estimators averaged over several runs. Number of runs N runs , number of particles 
N p , and time step At, have been chosen in such a way to balance time discretization 
error with statistical fluctuations. 

The first run consists to take a large number of particles and to make time av- 
eraging in order to minimize the fluctuation errors: iV mns = 10 3 , N p = 10 4 and 
At = 0.01. The total computational time to compute the evolution of moments is 
in this case 113 min. 23 sec. On the other hand we perform a second run where we 
use more averaging to minimize fluctuation errors: N rnns = 5. 10 3 , N p = 5. 10 3 and 
At = 0.01. The computational time is now 100 min. 

We remark that the computations have been obtained by using the hard spheres 
model (VHS with 7 = 1), which is the most realistic. In this case, the computational 
time of the Monte-Carlo methods becomes larger than the case when we consider 
pseudo-Maxwell molecules (VHS with 7 = 0) for which the collision kernel is con- 
stant and no rejection is needed. For the spectral method the computational cost is 
independent of the collision kernel. 
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Figure 5. 3D homogeneous case: time evolution of the kinetic en- 
tropy H and high-order moments Aii, M.q and Ais of f(t,v) for the 
fast and classical spectral methods, and the Monte- Carlo methods. 

From the comparison, it is obvious that, for three dimensional computations, the 
greater cost of the classical spectral scheme (with the same number of degrees of 
freedom) is compensated by a much greater accuracy, allowing better results with 
the same computational cost. Moreover, with the fast algorithm, the spectral scheme 
really becomes competitive in terms of computational time since the accuracy is not 
affected when we use a few number of discrete angles (for instance M = 4). In 
Figure we compare the accuracy on the evolution of high-order moments with 
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the different methods: the fourth order moments are very close, but the results 
obtained with the Monte-Carlo methods are affected by fluctuations on the tail of 
the distribution function, which are difficult to remove (see the evolution of the 8-th 
order moment). Note that for 2D and 3D pseudo-Maxwell molecules, comparisons 
had also been performed between Monte-Carlo methods and the classical spectral 
method in [HJ Section 6.3]. 

5.3. Stability of spectral methods with respect to non smooth data. In this 
subsection we perform some numerical simulations in order to study the behavior of 
the spectral methods when applied to non smooth data. We consider the following 
distribution 



and use a mollified initial datum, which suitably approximates moments. We per- 
form two numerical simulations using the fast spectral method with 32 2 and 64 2 grid 
points and plot the evolution of the entropy and the fourth order moment in Figure 
El Even for this discontinuous initial datum, we observe that for the two configu- 
rations the numerical entropy is decreasing and both numerical solutions converge 
to the same steady state. Moreover, we plot the evolution of the distribution func- 
tion with respect to time in Figure [7| the fast spectral method is very stable for 
this numerical test even if spurious oscillations are first generated, the distribution 
becomes smooth and converges to an approximated Maxwellian. As expected from 
a Fourier- Galerkin method, the accuracy degenerates in the discontinuity region. 
However, surprisingly the method seems to remain stable. 

Finally, let us mention that we have also performed some numerical tests when the 
initial datum approximates a Dirac distribution. The spectral method is still stable 
even if spurious oscillations are generated, but this problem is inherent to gridded 
methods and we refer to for a rescaling method, which follows the variation of 
the distribution function and allows to treat concentrated distributions. 



In |2U1 > we performed several numerical simulations to compare the spectral 
scheme with Monte-Carlo methods and showed that when we are interested in the 
transient regime, the deterministic method becomes very efficient. Obviously, the 
fast spectral method will still improve the computational cost. 

In the sequel, since we will be interested in the study of the trend to a global 
equilibrium state of the kinetic equation, we avoid the use of a splitting method 
by solving the whole non homogeneous equation in time by a second order Runge- 
Kutta method. Clearly the spectral methods apply straightforwardly to the collision 
operator also in this situation. The transport part is treated by the positive flux 
conservative method (see [11112111201 for further details). 




6. Application to the non homogeneous case 
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6.1. Definition of the problem. We consider the full Boltzmann equation in 
dimension d = 2 on the torus 

^ + v-V x f = Q(f,f), xe[0,L] 2 ,veR 2 

with periodic boundary conditions in x. We first introduce the hydrodynamical 
fields associated to a kinetic distribution f(t,x,v). These are the (d + 2) scalar 
fields of density p (scalar), mean velocity u (vector valued) and temperature T 
(scalar) defined by the formulas (|2.7j) . Whenever f(t,x,v) is a smooth solution 
to the Boltzmann equation with periodic boundary conditions, one has the global 
conservation laws for mass, momentum and energy 



d 
di 
d 
dt 



/(t, x, v) dxdv = 0, 

[0,L] 2 xR 2 

f(t, x, v) v dxdv = 0, 

[0,L] 2 xR 2 



d f \v\ 2 

— / f(t,x,v) dxdv — 0. 

dt J[0.Ll 2 xR 2 2 



'[0,L] 2 xR 2 

Therefore, without loss of generality we shall impose 

/ f(t,x,v)dxdv = l, / f(t,x,v)vdxdv 

J[0,L] 2 xIR 2 J[0,L] 2 xR 2 
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Figure 7. Time evolution of the distribution function for 64 x 64 grid points, 
and 

/ f(t.x.vV- 



12 



\V\ 

f(t, x, v) —— dxdv = 1. 



'[0,L] 2 xR 2 

These conservation laws are then enough to uniquely determine the stationary state 
of the Boltzmann equation: the normalized global Maxwellian distribution 

(6.1) Mg (v) = ^exp f M * 
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We shall use the following terminology: a velocity distribution of the form (J6.1j) will 
be called a Maxwellian distribution, whereas a distribution of the form 



will be called a local Maxwellian distribution (in the sense that the constants p, u and 
T appearing there depend on the position x). We also define the notion of relative 
local entropy Hi, the entropy relative to the local Maxwellian, and the relative global 
entropy H g , the entropy relative to the global Maxwellian distribution, by 



Our goal here is to investigate numerically the long-time behavior of the solution 
/. If / is any reasonable solution of the Boltzmann equation, satisfying certain 
a priori bounds of compactness (in particular, ensuring that no kinetic energy is 
allowed to leak at large velocities), then it is possible to prove that / does indeed 
converge to the global Maxwellian distribution M g as t goes to +00. Of course, 
obtaining these a priori bounds is extremely difficult; as a matter of fact, they have 
been established only in the spatially homogeneous situation (which means that the 
distribution function does not depend on the position variable x, see the survey in 
[35] ) or in a close-to equilibrium setting (see in particular j2H] for the torus), and it 
still constitutes a famous open problem for spatially inhomogeneous initial data far 
from equilibrium. More recently, Desvillettes and Villani JB], Guo and Strain [37] 
were interested in the study of rates of convergence for the full Boltzmann equation. 
Roughly speaking in [IB] , the authors proved that if the solution to the Boltzmann 
equation is smooth enough and satisfies bounds from below of the form 



Vt > 0, x G [0, L} 2 , v G M 2 , f(t, x, v) > K e~ A ° H?0 (A , K > 0, q > 2), 



(although this bound can be shown to be a consequence of the regularity bounds, 
see [20]) then (with constructive bounds) 



which means that the solution converges almost exponentially fast to the global 
equilibrium (namely with polynomial rate 0(t~ r ) with r as large as wanted). 

The solution / to the Boltzmann equation satisfies the formula of additivity of 
the entropy: the entropy can be decomposed into the sum of a purely hydrodynamic 
part, and (by contrast) of a purely kinetic part. In terms of H functional: one can 
write 



(6.2) 





\\f(t) - M g \\ = o(r~) 
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In fact, we can also show that 

Hi(t) < H g (t), Vi > 0. 

Moreover, the Csiszar-Kullback-Pinsker inequality asserts that (when the total mass 
of the solution is normalized to 1) 

H{f\M)>\\\f-Mf Li . 

In other words, controlling the speed of convergence of the entropy to its equilibrium 
value is enough to control the speed of convergence of the solution to equilibrium, 
in very strong sense. 

Moreover in jTHj, Desvillettes and Villani conjectured that time oscillations should 
occur on the evolution of the relative local entropy. In fact their proof does not rule 
out the possibility that the entropy production undergoes important oscillations in 
time, and actually most of the technical work is caused by this possibility. 

6.2. Description and interpretation of the results. Here, we performed sim- 
ulations on the full Boltzmann equation in a simplified geometry (one dimension 
of space, two dimensions of velocity, periodic boundary conditions, fixed Knudsen 
number) with the fast spectral method to observe the evolution of the entropy and 
to check numerically if such oscillations occur. Clearly this test is challenging for 
a numerical method due to the high accuracy required to capture such oscillating 
behavior. 

Then, we consider an initial datum as a perturbation of the global equilibrium 

M g 

(6.3) fo(x,v) = — (1 + A cos(k x)) exp(— |t>| 2 /2), x G [Q,L],v G R 2 
In 

for some constants A > and k = 2ir/L. 

In Figures |H1 and 03 we are indeed able to observe oscillations in the entropy pro- 
duction and in the hydrodynamic entropy. The strength of the oscillations depends 
a lot on the length L of the domain, which is consistent with the fact that such 
oscillations are never observed in the spatially homogeneous case (L = 0). The 
superimposed curves yield the time evolution respectively of the total H functional 
and of its kinetic part. In all local Maxwellian distribution is chosen for 

initial datum; the first plot corresponds to L — 1 and the second one to L — 4. 
Some slight oscillations can be seen in the case L — 1, but what is most striking is 
that after a short while, the kinetic entropy is very close to the total entropy: an 
indication that the solution evolves basically in a spatially homogeneous way (con- 
trary to the intuition of the hydrodynamic regime). On the contrary, in the case 
L = 4, the oscillations are much more important in frequency and amplitude (note 
that this is a logarithmic plot): the solution "hesitates" between states where it is 
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very close to hydro dynamic, and states where it is not at all. Further note that the 
equilibration is much more rapid when the box is small, and that the convergence 
seems to be exponential. 

It is in fact possible to give a simple interpretation of these oscillations thanks to 
the work jTTj . Since this effect is observed near the global equilibrium one can replace 
the Boltzmann collision operator by the linearized Boltzmann collision operator 
(moreover the oscillation effect is effectively observed for the linearized Boltzmann 
collision operator as well). Then it is straightforward that the computations above 
correspond to observing the time evolution of one Fourier mode in x (here with 
frequency k and amplitude A ). Hence by an obvious rescaling, this evolution is 
given by the semi-group T ko / L (t), where T k {t) is defined in [T7] (this is the semi-group 
for the k-th Fourier mode in x for the linearized equation). An asymptotic study 
of the spectrum of its infinitesimal generator for small frequencies k was done in 
[T7j . The dominant term in terms of long time behavior (i.e., the one with the lower 
rate of decrease) is given by the (d + 2) "hydrodynamical eigenvalues" . Moreover 
explicit computations are available for the expansions of these eigenvalues according 
to e = \k\ near k = 0. 

At first order in e, the eigenvalues vanish, except for two of them, which are purely 
imaginary. They are given by 

' h=ie^l + 2/d + 0(e 2 ), 
< h = -te^l + 2/d + 0(s 2 ) 1 
K h = ■ ■ ■ = Id+2 = 0. 

Therefore for \k \/L << 1 (realized for instance when fc is fixed and L is large 
enough), this analysis gives us the dominant imaginary term in the eigenvalues. In 
this regime, one should thus observe oscillations with frequency ^/l + 2/d \k \/L. 
Thus the period of oscillations should be given by 27i(l + 2/d)-^ 2 L/\k \, which can 
be checked with the numerical simulations. Indeed, in Tabled we give the ratio of 
the period of oscillations uo with the length box L. The numerical results agree well 
with the analytical computations u/L~ l/v2. 

We also observe that the damping rate is related to the length box and is pro- 
portional to 1/L 2 when L becomes large (see Tabled a L 2 ~ constant). This is 
coherent with the fact that no real value occurs in the "hydrodynamical" eigenvalues 
until the second order in e = \k\. The coefficients for the order 2 in the expansion 
are computed in [T7] ; they are purely real and they can be expressed simply in terms 
of the dimension d, the viscosity coefficient n and the heat conductivity A of the gas 
(indeed these coefficients are related with the Navier-Stokes limit of the Boltzmann 
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equation). Namely they are given by 

i?l = i?2 = 
R 3 = ■ ■ • = 

Rd+2 



A 



d + 2 

Rd+l = 

Xd 



V 
2 



r\ d 



2(d-l) 



d + 2 



Therefore for \ko\/L << 1 , the damping rate is given by the minimum among these 
values. 

To conclude these tests, we performed a last numerical experiment to evaluate 
the robustness of the theory of the trend to equilibrium. We have chosen an initial 
datum which is far from the equilibrium 
(6.4) 

1 



fo(x,v) 



[I + A cos(k x)) 



exp 


( 


v - v 


2 \ 


+ exp 


< 


v + v 


2 \ 


K Hh ) 


\ z ' tJ th J 



with d = (1/2, 1/2) and v t h = v3/2. We present the time evolution of the relative 
entropies Hi and H g in log scale and observe that initially the entropy is strongly 
decreasing and when the distribution function becomes close to a local equilibrium, 
some oscillations appear with the good frequency uj/L = l/y/2 and damping rate a 
(see Figure EI)- 

Remarks: 

1. Now numerical methods for the Boltzmann equation -such as the one presented 
in this paper- become able to provide very accurate simulations of the transient 
regime towards equilibrium with reasonable cost, even in the inhomogeneous case. 
This could be used to explore numerically the spectrum of the linearized Boltzmann 
collision operator (in the homogeneous case), or, more interestingly, the spectrum 
of the linearized Boltzmann collision operator together with the transport term in 
the inhomogeneous case. For instance the exponential rate of convergence is directly 
readable on the figures above, and provides a numerical estimation of the spectral gap 
(that is the real part of the first non-zero eigenvalue) of this operator (it is known 
since Ukai jlH] that this operator has a spectral gap in the torus, see also [T3]). 
Moreover by a frequency analysis of the curve of the time evolution of the relative 
entropy or the L 1 distance to the equilibrium, it could be possible also to describe 
other eigenvalues: as long as they have different imaginary part, it should be possible 
(in principle) to extract from the frequency analysis the curve corresponding to their 
contribution in the evolution semigroup, and then to compute their real part which 
corresponds to the exponential rate of decay of this curve. 
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Length box 


oscillation frequency uj 


uj/L 


damping rate a 


-aL' 2 


L = ir/2 


01.10 


0.701 


-6.521 


16.04 


L= 7T 


02.25 


0.716 


-2.202 


21.71 


L = 2vr 


04.50 


0.716 


-0.641 


25.26 


L = 3tt 


06.61 


0.701 


-0.285 


25.31 


L = 4vr 


08.78 


0.699 


-0.160 


25.27 


L = 8tt 


17.57 


0.699 


-0.040 


25.35 



Table 4. Influence of the length box: damping rate and oscillation 
frequency for the relative entropy with respect to the local Maxwellian 
Hi{t) using 64 x 64 x 64 with A = 0.1 in Q . 



2. A recent work j2Z| gave a detailed pointwise study of the Green function 
for the linearized Boltzmann equation in the domain x G Q = M. In particular 
in this case, there study shows that the long-time behavior is governed by "fluid- 
like waves" (corresponding to the waves of the linearized Euler and Navier-Stokes 
equations) whose amplitude decreases polynomially, whereas the amplitude of the 
"kinetic part" of the Green function decreases exponentially. We think it likely that 
this study could be extended to the torus, where the amplitude of the fluid and 
kinetic parts of the Green function should both decrease exponentially. Moreover 
the rate of decay of the kinetic part should not depend on the size of the box, 
whereas the rate of decay of the fluid part should do. Hence for a box small enough, 
the long-time behavior should be governed by the kinetic part of the Green function 
(that is like the spatially homogeneous Boltzmann equation), whereas for a box big 
enough, the long-time behavior should be governed by the "fluid-like waves" . This is 
precisely what we observe numerically, and thus this theoretical study could provide 
a rigorous proof of the numerical observations above, at least in the linearized regime. 

7. Conclusions 

In this paper we have introduced and deeply tested a class of new fast algorithms 
for the computation of the Boltzmann collision operator. These methods allow to 
reduce the computational cost from 0(n 2 ) to 0(nlog 2 n). We give computational 
evidence of the great performance of the schemes which can provide a dramatic 
speed up in computing time of deterministic schemes by making them competitive 
with Monte-Carlo methods where higher accuracy is required. A first numerical 
application to a non trivial problem in the space non homogeneous case confirms 
the strong computing potential of the new schemes. 

Other methods such as singular value decomposition, fast multipole methods [9 , 
separated representations in high- dimensional problems (see works by G.Beylkin for 
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Figure 8. Influence of the length box: relative entropy with respect 
to the local Maxwellian Hi(t) using 64 x 64 x 64 for L= n; 2tt; 3tt and 
4tt with A = 0.1 in (IQl . 



instance P), wavelets, . . . could have been used to search for some decomposition 
of the form (j4.19|) . However to our knowledge it is not known at now how to obtain 
the properties described above on this decomposition with these methods. 
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